subroutine  compute_Richter_ML(c,nd,stanum,evnum,nsta,nev,sigma2_RH)
    implicit none
    
    include 'seisev_EVENT_PAHSE_structure.FI'
    include 'inversion_of_residuals_ML.FI'
    
    integer i,k,nd,isw,nv
    integer stanum(200000),evnum(200000),nsta,nev
    integer RichDist(100),first,last,dst
    real*4  RichVal(100),dif,corr
    real*4  v(5000),wmag,wsd,wmax,wmin,ww(5000),ruse
    real*8  sigma2_RH
    data isw/0/

   
   ! compute (not read) epicentral distances from hypocentral and depth.
    do i=1,nd
       c(i).MLRH_dist = SQRT(  c(i).s_h_hypdist * c(i).s_h_hypdist -  &
                               c(i).depth       * c(i).depth )
    enddo


   ! compute all the station Richter magnitudes
    if(isw .eq. 0)then
       isw = 1
       call read_Richter_tab(nv,RichDist,RichVal)
    endif
    
    do i=1,nd       
       dst = NINT(c(i).MLRH_dist)
       ! compute correction
       call compute_Richter_corr(dst,corr,nv,RichDist,RichVal)
       ! compute station magnitude
       c(i).MLRH_sta = log10(c(i).WA_pha_amp) + corr
    enddo
    
    
    ! go throught the events and compute the event Richter Magnitudes.
    

    first = 1
    last  = 1
    do while( evnum(last) .eq. evnum(first) )
      last = last +1
    enddo
    last = last -1

    sigma2_RH=0.0d00
    do while (last .le. nd  .and. first .le. nd)

      k= 0
      do i=first,last
          k=k+1 
          v(k) = c(i).MLRH_sta
      enddo

      ruse = 0.3
      wmag = 0.
      wsd = 0.
      call Huber_dev(k,v,ruse,wmag,wsd,wmax,wmin,ww)
      ! memorizza i risultati   !
      !
      

      k=0
      do i=first,last
          k=k+1
          c(i).MLRH_ev = wmag     !event mag
          dif = c(i).MLRH_sta - wmag
          sigma2_RH = sigma2_RH + dif*dif
      enddo   
                 
      ! vai al prossimo evento
      first = last + 1
      last = first
      do while( last .le. nd .and.  evnum(last) .eq. evnum(first))
        last = last +1
      enddo
      last = last -1

    enddo        
    
    write(*,*)' calcolo delle ML Richter. Sigma2 = ',sigma2_RH    
    pause


return
end
    



subroutine compute_Richter_corr(dst,corr,nv,RichDist,RichVal)
! ATTENZIONE: questa routine funziona solo se la tabella
! inizia dalla distanza 0.

integer dst
real*4  corr,dd
integer nv,i,RichDist(100)
real*4       RichVal(100)

dd=float(dst)

do i=1,nv

  if(RichDist(i) .gt. dd)then
    a = RichDist(i) - RichDist(i-1)
    b = RichVal(i) - RichVal(i-1)
    corr = RichVal(i-1) + b/a*(dd - RichDist(i-1))
    return
  endif 
enddo
! oltre la distanza massima:
corr = RichVal(nv)
return
end
    